Bond-order-wave phase and quantum phase transitions 
in the one-dimensional extended Hubbard model 

Pinaki Sengupta 

Department of Physics, University of Illinois at Urbana- Champaign, 1110 West Green Street, Urbana, Illinois 61801 

Anders W. Sandvik 

Department of Physics, Abo Akademi University, Porthansgatan 3, FIN-20500, Turku, Finland 

David K. Campbell 

Departments of Physics and of Electrical and Computer Engineering, Boston University, 44 Cummington Street, 

Boston, Massachusetts 02215 
(February 1, 2008) 

We use a stochastic series expansion quantum Monte Carlo method to study the phase diagram of 
the one- dimensional extended Hubbard model at half filling for small to intermediate values of the 
on-site {U) and nearest-neighbor {V) repulsions. We confirm the existence of a novel, long-range- 
ordered bond-order-wave (BOW) phase recently predicted by Nakamura (J. Phys. Soc. Jpn. 68, 
3123 (1999)) in a small region of the parameter space between the familiar charge-density-wave 
(CDW) state for 1/ > C//2 and the state with dominant spin-density-wave (SDW) fluctuations 
for V < U/2. We discuss the nature of the transitions among these states and evaluate some 
of the critical exponents. Further, we determine accurately the position of the multi-critical point, 
{Um, Vm) = (4.7±0.1, 2.51±0.04) (in energy units where the hopping integral is normalized to unity), 
above which the two continuous SDW-BOW-CDW transitions are replaced by one discontinuous 
(first-order) direct SDW-CDW transition. We also discuss the evolution of the CDW and BOW 
states upon hole doping. We find that in both cases the ground state is a Luther-Emery liquid, i.e., 
the spin gap remains but the charge gap existing at half-flUing is immediately closed upon doping. 
The charge and bond-order correlations decay with distance r as r~^'' , where Kp is approximately 
0.5 for the parameters we have considered. We also discuss advantages of using parallel tempering 
(or exchange Monte Carlo) — an extended ensemble method that we here combine with quantum 
Monte Carlo — in studies of quantum phase transitions. 
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I. INTRODUCTION 

The one-dimensional (ID) extended Hubbard model 
has been extensively studied in recent years, both as an 
important theoretical test-bed for studying novel con- 
cepts in ID (e.g., spin-charge separation) and methods, 
(e.g., quantum Monte Carlo, exact diagonalization, and 
the density matrix renormalization group) and as a use- 
ful model for several classes of quasi ID materials in- 
cluding copper-oxide materials related to thephigh-Tc 
cuprate superconductors Jil conducting polymersB and or- 
ganic charge-transfer salts.cl General ID extended Hub- 
bard models differ from the standard Hubbard model, 
which includes only an on-site electron-electron inter- 
action ([/), by the addition of longer-range interactions 
which are necessary to explain several experimentally ob- 
served effects in real materials, e.g., excitons in conduct- 
ing polymers. The simplest extended Hubbard model 
(henceforth, EHM), on which we focus in this article, 
consists of adding a nearest neighbor interaction V . If 
the interaction parameters are assumed to arise solely 
from Coulomb interactions, both U and V are repulsive 
(positive), and U > V . However, viewed as phenomeno- 
logical parameters incorporating the effects of additional 



(e.g. electron-phonon) interactions, the ranges of these 
parameters can be much broader, including U,V < 0. 
The Hamiltonian is 

H = -tJ2i4+l.aC^.a + h.C.) 

+c^EKT-5)Ki-^) 

i 

-fyE(n,+i - l)(n. - 1) -M^n,, (1) 

i i 

where c] ^{ci^a) creates (annihilates) an electron with 
spin a at site i, t is the hopping integral between ad- 
jacent sites and fi is the chemical potential. Henceforth 
we set t — 1 and express the interaction parameters U 
and V in units of t. 

The ground state phase diagram of the EHM at half 
filling (/i — 0) has been extensively studied using both 
analytical and numerical methods. Despite the appar- 
ent simplicity of the model, the phase diagram shows 
surprisingly rich structure. In the limit V—0 (the stan- 
dard Hubbard model), the Hamiltonian (|l|) can be dk 
agonalized exactly using the generalized Bethe Ansatz.cl 
For V ^ 0, the model has been studiecLusiag perturba- 
tive methods and numerical simulations .Erila Broadly, the 
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phase diagram consists of insulating phases with domi- 
nant charge-density-wave (CDW) and spin-density-wave 
(SDW) characters and metalHc phases where singlet and 
triplet superconducting correlations dominate. In the 
physically relevant region for "Coulomb-only" parame- 
ters {U, V > 0), the system is in a CDW phase for large 
V/U and in an state with dominant SDW fluctuations for 
small V/U. The CDW phase has broken discrete sym- 
metry characterized predominantly by alternating dou- 
bly occupied and empty sites and exhibits long-range 
order. The SDW phase, on the other hand, has con- 
tinuous symmetry and hence cannot exhibit long-range 
order in ID (by the Mermin- Wagner theorem). Instead, 
it is a critical state characterized by the slow (algebraic) 
decay of the staggered spin-spin correlation function. In- 
deed, in the limit U ^ 1, U ^ V, the model reduces 
to an effective Heisenberg model with J ^ 1/{U — V). 
For small U and V {U,V <C 1), the boundary between 
the CDW and the SDW phases was predicted to be at 
U — 2V using weakxgupling renormalization group tech- 
niques (" g-ology" ) .Q'LI Strong coupling calculations us- 
ing second-order perturbation theory also gave the same 
phase boundary {U = 2V) between the CDW-and the 
SDW phases for large U and V {U,V :^ l)m For in- 
termediate values of the parameters, the phase boundary 
was found to be shifted slightly away from the U = 2V 
line such that the SDW phase is enlaanced, as shown 
by quantum Monte Carlo simulationaj'Q as well as strong 
coupling calCjUlations using perturbation theory up to the 
fourth order .Ej Moreover, the nature of the transition is 
quite different in the two coupling regions, changing from 
continuous (second-order) in the weak coupling limit to 
discontinuous (first-order) in the strong coupling limit. 
Estimates for the location of the multi-critical point, 
where the nature of the transition changes, have pa nged. 
from U,n ^ 1.5 to U,n ^ 5 (and K„ « f7™/2)j3"By 
Despite the broad uncertainty in the actual value of the 
tricritical point, the phase diagram was believed to be 
well understood. 

Recently, however, by studying the EHM ground state 
broken symmetries using level crossings in excitatioa 
spectra obtained by exact diagonalization, NakamuraE^ 
has argued for the existence of a novel bond-order-wave 
(BOW) phase for small to intermediate values of U and V 
in a narrow strip between the CDW and the SDW phases. 
The BOW phase is characterized by alternating strengths 
of the expectation value of the kinetic energy operator on 
the bonds. It is predicted to be a state where the dis- 
crete (two-fold) symmetry is broken and should hence ex- 
hibit true long-range order. Nakamura thus argues that 
the transition between CDW and SDW phases in this 
region is replaced two separate transitions: (i) a continu- 
ous transition from CDW to BOW; and {ii) a Kosterlitz- 
Thouless spin-gap transition from BOW to SDW. The 
BOW region vanishes at the multi-critical point beyond 
which the transition between CDW and SDW phases is 
direct and discontinuous. A schematic phase diagram 
including Nakamura's BOW state is shown in Fig. m 



V 




FIG. 1. Schematic ground state phase diagram of the 
EHM at half fiUing, as proposed by Nakamura. The CDW 
and BOW phases are long-range-ordered (broken-symmetry) , 
whereas the SDW phase has no broken symmetry but exhibits 
an algebraically decaying spin-spin correlation function. 

Considering the long history of the ID EHM and the 
large number of studies of the U « 2V region with a va- 
riety of analytical and numerical tools, the proposal of 
a new phase is certainly remarkable. Importantly, the 
level crossing method used by Nakamura cannot by itself 
exclude the conventional scenario of a direct SDW-CDW 
transition for the whole range oi U,V > 0; a level jCkoss- 
ing corresponding to this transition was also foundliS be- 
tween the SDW-BOW and BOW-CDW crossing curves. 
The position of the BOW-CDW level crossing is, how- 
ever, in closer agreement with the strong-coupling result 
for the vanishing of the CDW order, and this was taken 
as evidence of a long-range-ordered BOW in the ground 
state for certain parameters. It is important to confirm 
this hitherto undiscovered phase using other methods. 

To attempt this confirmation, we have used the 
highly efficient stochastic fsjip^ expansion (SSE) quan- 
tum Monte Carlo methodE-lil3 to study the EHM at 
half filling in the vicinity oi U = 2V. This method al- 
lows us to probe directly the spin- charge- and bond- 
order correlations in the ground state of lattices with 
more than one hundred sites (up to 256 sites were used 
in this study). Using finite-size scaling techniques for 
the various order parameters, we confirm the existence 
of a BOW state with spin and charge gaps in a region 
very close to that predicted by Nakamura for small U, V. 
We also further improved the SSE simulations by ap- 
plying a quantum version of the thermaLj^atallel tem- 
pering scheme (or exchange Monte Carlo )tJL2l for simu- 
lations close to and across the phase boundaries. This 
"quantum parallel tempering" greatly reduced the ef- 
fects of "sticking" — where the simulation gets trapped 
in the wrong phase close to a phase boundary — and 
was found to be particularly useful for the discontinuous 
(first-order) direct SDW-CDW transition. As a conse- 



2 



quence, we were able to obtain a more accurate estimate 
for the location of the multi-critical point Vm) where 
the BOW phase vanishes and is replaced by a first-order 
SDW-CDW transition line. As we discuss below, we find 
Um = 4.7 ± 0.1, V„, = 2.51 ± 0.04. 

In order to investigate the possibility of soliton lattices 
forming out of the long-range CDW and BOW states 
when doping away from half-filling, we have also carried 
out some simulations of lightly doped systems. We find 
that the in both cases the ground state is a Luther-Emery 
liquid, with a spin gap and slow algebraic decay (~ r~^p , 
with Kp « 0.5) of the dominant CDW and BOW corre- 
lations. 

The remainder of the paper is organized into four sec- 
tions and two appendices. In Sec. |l|we briefly sketch the 
SSE method and introduce the different observables we 



study. In Sec. |I| we present the results of our simulations 
at half-filling and discuss their interpretation. Doped sys- 
tems are considered in Sec. IV. We conclude with a brief 
summary in Sec. In Appendix A we present some 
important details of the extension of the SSE method 
to allow efficient loop updates for fermions. We illus- 
trate the advantages of the quantum parallel tempering 
scheme in Appendix B. 



with "trapping" close to a first-order phase transition, 
i.e., the simulation can get stuck in the wrong phase very 
close to the critical point. There are also problems with 
slow dynamics in long-range ordered phases with a bro- 
ken discrete symmetry (such as BOW or CDW phases). 
In order to overcome these problems we have developed 
a "quantum parallel tempering" scheme — a geneijaliza- 



tion of the thermal parallel tempering methoccj til com- 
monly used to equilibrate classical spin glass simulations. 
The method amounts to running several simulations on a 
parallel computer, using a fixed value of U and different 
closely spaced values of V at and around the critical value 
Vc- Along with the usual Monte Carlo updates, we at- 
tempt to swap the configurations for processes with adja- 
cent values of V at regular intervals (typically after every 
Monte Carlo step) according to a scheme that maintains 
detailed balance in the space of the parallel simulations, 
as explained in Appendix B. In contrast to Ref. we 
here find parallel tempering to be particularly useful in 
the study of the first-order transition, where the prob- 
lem of trapping is the most pronounced. In Appendix B 
we also present a comparative example to illustrate the 
improvement obtained by parallel tempering. 



II. NUMERICAL METHODS AND 
OBSERVABLES 

A. The SSE method and its fermion loop-update 
extension 

The SSE method0'0 is a finite-temperature quantum 
Monte Carlo method based on importance sampling of 
the diagonal elements of the Taylor expansion of e~^^ , 
where (3 is the inverse temperature /3 = t/T. Ground 
state expectation values can be obtained using suffi- 
ciently large values of /3, and there are no approximations 
beyond statistical errors. Recently, in the context of spin 
systemsjlj an efficient "operator loop update" was de- 
veloped to sample the operator sequences appearing in 
the expansion. The resulting method has_pp3ven to be 
very efficient for several different models.t3t3 To apply 
the most efficient variant of SSE method to the EHM, 
we need to generalize the previous operator loop update 
scheme to spinful fermions. This is an important exten- 
sion, but because of its technical nature we have relegated 
our detailed discussion of it to an appendix. 

We have applied the SSE algorithm to the ID EHM 
for system sizes ranging from = 8 to 256 sites, with 
maximum inverse temperatures (3 chosen appropriately 
to isolate the ground state. We have verified the correct- 
ness of the simulation code by comparing A^ = 8 results 
with exact diagonalization (Lanczos) results. 

Although the operator-loop update is indeed signifi- 
cantly more efficient than previous local updates for sam- 
pling of the SSE configurations, we still have problems 



B. Observables 

In addition to the ground state energy, E = {H)/N, 
the observables we study include the static structure fac- 
tors and susceptibilities corresponding to the different 
phases (CDW, SDW and BOW). The structure factors 
are given by 



5sw(9) = ^E^''^'"'^('^J^fc)' 

ScDwiq) = ^ - K)': 

SBOwiq) = ^ E ^'"^'^'^ (^^■^'») - (^^■)'' (2) 

j,k 



where 



= E -t- h.c.) 



(3) 



is the kinetic energy operator associated with the jth 
bond. The corresponding static susceptibilities are given 

by 

XSDwiq) = ^E^*'''"'^ f dr{S^^{r)Sm) (4) 

]:k 

and analogous expressions for xcDw{<l) and xbow{q)- 
Since all the phases mentioned have a period 2, the stag- 
gered structure factor and susceptibilities are the most 
important observables. We define order parameters for 
the phases in terms of the staggered structure factors: 
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where a — CDW, SDW, or BOW. We have also studied 
the charge stiffness constant, pc- It is defined as the 
second derivative of tlie internal energy per site, E, with 
respect to a twist, (pf^ 



(6) 



under which the hopping term in the Hamiltonian (|l|) is 
replaced by 



(7) 



The spin stiffness constant, ps, is defined by a similar 
expression, with the hopping term now being replaced 
by 

= -i^(e-^^"c]^_i,,c,, . + h.c), (8) 

with 01 = —(f>i — <j). In the framework of the SSE 
method, the estimators for the charge and spin stiffness 
are given in terms of expectation values of squared wind- 
ing numbers (see Appendix A). 



III. RESULTS AT HALF-FILLING 

As noted above, we have studied chains-with N up 
to 256 with periodic boundary conditions E3 Typically, 
an inverse temperature of = 2N was sufficient for the 
calculated properties to have converged to their ground 
state values, except in the case of iV = 256, for which 
(3 = AN was needed for some quantities. In this section 
we first discuss our evidence for the existence of a long- 
range BOW phase, then our analysis of the continuous 
BOW-CDW and SDW-BOW transitions for smaU [U, V), 
the discontinuous SDW-CDW transition for large (J7, V"), 
and finally our determination of the location of the multi- 
critical point separating these transitions. 



A. Existence of the BOW phase 

Plots of the variation of the staggered susceptibilities 
corresponding to the three different phases — CDW, 
SDW, and BOW — show the existence of strong BOW 
fiuctuations in a region with V ~ 11/2 in parameter space 
where Nakamura predicted a BOW state. Fig. |^ is such 
a plot for ?7 = 4 and 1.7 < V < 2.3. In a long-range 
ordered phase (BOW, CDW), the corresponding xi''^) is 
expected to diverge with increasing system whereas the 
other two susceptibilities should converge to constants. 
In the SDW phase there is no long-range order but alge- 
braically decaying correlations of both SDW and BOW 




FIG. 2. The variation with V (at fixed (7 = 4) of the 
staggered susceptibilities (CDW, BOW, and SDW, from the 
top) in the neighborhood of the BOW phase predicted by 
Nakamura (the vertical dashed lines show the predicted 
SDW-BOW and BOW-CDW boundaries). The statistical 
errors are typically of the order of the size of the symbols 
(slightly larger for the N = 128 CDW at high V). The scans 
for TV = 16 and 32 were obtained in single parallel tempering 
simulations, whereas those for = 64 and 128 consisted of 
two and four non-overlapping runs, respectively. 



nature; hence xsDwiT^) and XBOwiT^) should both di- 
verge here, but the BOW divergence should be much 
slower than in the long-ranged BOW phase. These be- 
haviors are indeed seen in Fig. ^ with the susceptibilities 
for SDW, BOW, and CDW dominating in turn as V is in- 
creased. The BOW-CDW phase boundary can be quite 
well resolved since it involves a standard second order 
(continuous) phase transition. On the other hand, the 
SDW-BOW boundary is more difficult to locate, for it in- 
volves a Kosterlitz-Thouless transition in which the spin 
gap opeps exponentially slowly as one enters the BOW 
phasejl3 resulting in only a slow decay of the staggered 
SDW susceptibility in the BOW phase for the system 
sizes accessible in our work. 

Fig. H shows ln[xa(7r)] and ln[S'Q(7r)] versus ln[A^] for 
the parameters [U^V] = (4,2.14) for which the ground 
state should be inside the BOW phase. We find that 
both XBOwiT^) and SBOwi^^) diverge strongly with sys- 
tem size, whereas the structure factor and susceptibility 
corresponding to CDW have a maximum and then de- 
crease with system size for large N . The SDW structure 
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FIG. 3. ln[x(7r)] and ln[S'(7r)] vs ln[iV] for the diflFerent 
phases at (7 = 4, V — 2.14 and system sizes A*' up to 256. 
The dashed line in the ^(Tr) panel has slope 1. 

factor appears to have converged for TV = 256 but the sus- 
ceptibiUty still shows a weak growth — in a spin-gapped 
BOW phase it should eventually converge, too, but if the 
gap is very small the convergence occurs only for much 
larger systems. The growth with N seen here is much 
slower than N, which should be the asymptotic behav- 
ior in an SDW phase for any spin-rotationally invariant 
ID systemj23 and the growth slows with increasing N. 
Hence an asymptotic divergence of xsDwi''^) can be ex- 
cluded. The dominant asymptotic characteristic of the 
ground state is clearly BOW. The system sizes consid- 
ered are not large enough for SBowi''^) to have reached 
the asymptotic behavior ^ N expected if there is long- 
range order, which we will explain further below. The 
very fast divergence of xbow (tt) is expected on account 
of the two-fold degenerate BOW ground state. For finite 
N this degeneracy is not perfect, but an exponentially 
fast closing of the gap between the symmetric and an- 
tisymmetric linear combinations of the two asymptoti- 
cally degenerate symmetry-broken ordered states can be 
expected, which would eventually cause XBOwi''^) to di- 
verge exponentially. 

The most direct evidence for a long-range BOW comes 
from the the real-space kinetic energy correlation func- 
tion 
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As seen in Figure ^, this correlation function oscillates 



FIG. 4. The real-space BOW correlation function at (7 = 4, 
V = 2.14 for system sizes iV = 128 and 256. 



with period 2 and its magnitude decays considerably for 
short distances. For long distances there is a conver- 
gence to a constant, non-zero magnitude, which is the 
same within statistical errors for N = 128 and 256. The 
significant enhancement of the correlations at short dis- 
tances explain the deviations from the expected asymp- 
totic linear scaling of the integrated correlation function, 
Sbow{t^), for the system sizes shown in Figure ^. 

Further proof of the existence of the BOW phase is 
obtained by looking for spin and charge gaps in this 
region. Instead of calculating the gaps directly, which 
can not easily be done to high accuracy for large sys- 
tem sipes, we use the following indirect method: It is 
knowita that if the ground state of a ID system is gap- 
less in the spin sector, the Luttinger liquid parameter 
K„ governing the asftopptotic equal-time spin corijelation 
function is = l.liJ It has been further shownEj that 
the slope SsdwW) 1 1 gives K^j-K in the limit g — > 0. 
Hence, SsDw{q)l<l — > l/i" as g ^ 0. On the other hand, 
if the ground state has a spin gap, SsdwW) I <! ^ as 
q — > 0. With this criterion even a very small spin gap 
can be detected, since it is in practice sufficient to see 
that TrSsDw{q)/q decays below 1 for small q to con- 
clude that Kcr ^ 1 and hence that a spin gap must be 
present. Similarly, for a ground state with no charge 
gap, T^ScDwil) 1 1 Kp as g — > whereas if the ground 
state does have a charge gap, ScDw{q) / q — > as g ^ 0. 
Unlike where the value is fixed at 1 for spin rota- 
tionally invariant systems, the Luttinger liquid charge 
correlation parameter Kp is a function of U and V , and 
its precise_value for given U and V is not known (except 
a.t V = 0E3). Due to the logarithmic corrections typical 
for ID systems, it is very difficult to obaepie numerically 
that S s DW (q) / Q becomes exactly LLiTu Empirically, 
we have found that in the gapless case the value 1 is al- 
ways approaihed from above (which is the case also for 
spin systemsE3), and hence the detection of the spin gap 
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FIG. 5. SsDw{<l)/q and ScDw{q)/q vs q ior U = 4 and 
V = 2.14 and V = 1.8 {N = 128). 

using this quantity is not hampered by the log-corrections 
— if T^SsDw{l) / 1 decays below 1 one can conclude that 
here is a gap. 

Fig. H shows T^SsDwil)/ 1 and ttS'cdvkI'?)/'? versus q/ir 
for U — A and two values of V . One of the points 
iy = 2.14) is inside the BOW phase, whereas the other 
{V ~ 1.8) is in the SDW phase. The T^SsDw{<l)/q curve 
for V — 1.8 is close to 1 for a wide range of q values, 
whereas the V — 2.14 curve exhibits a sharp drop as 
g ^ indicating, respectively, the absence and the pres- 
ence of a spin gap. Similarly, the evidence for a vanishing 
limit of ScDW (9) and hence of a charge gap for V — l.% 
is clear. Since the point V = 2.14 is quite close to the 
critical point {Vc = 2.16) where the charge gap vanishes, 
the magnitude of the gap is very small and we need to 
go to still smaller g, i.e., larger system size, to see a pro- 
nounced effect like that for V = 1.8. Nevertheless, the 
downturn for the smallest g is a good indication of a gap. 

The opening of spin and charge gaps can also be de- 
tected in the spin and charge stiffness constants, which 
should vanish as A'^ — » 00 if there are gaps. The asymp- 
totic charge stiffness should hence be non-zero only ex- 
actly at the BOW-CDW phase boundary. The spin stiff- 
ness should be non-zero in the SDW phase, should ap- 
proach a constant value exactly atr4he phase boundary 
(with logarithmic size-corrections), Eil and vanish inside 
the CDW phase. In Fig. ^ we show the stiffness constants 
for J7 = 4 in the neighborhood of the BOW phase. As 
expected, the charge stiffness peaks at the BOW-CDW 




FIG. 6. Behavior of the charge and spin stiffness across the 
BOW-CDW boundary for U = A. The upper(lower) panel 
shows the charge (spin) stiffness. The vertical dashed lines 
indicate the position of the phase boundaries according to 
Nakamura. 

phase boundary and decreases rapidly away from it, con- 
firming the vanishing of the charge gap only at the phase 
boundary. The peak becomes very sharp for large sys- 
tem sizes, and the finite-size corrections to its location are 
small. We find this the most accurate way to locate the 
BOW-CDW phase boundary. The spin stiffness is clearly 
zero in the CDW phase, and a sharp decrease with in- 
creasing N is also seen for V values well inside the BOW 
phase. Since the spin gap opens up exponentially slowly 
at the SDW-BOW boundary it is difficult to locate the 
transition this way. Our data nevertheless indicate that 
the BOW phase at t/ = 4 may not extend down to the 
value V K, 1.82 obtained by Nakamura. We will discuss 
this phase transition and determine the transition point 
more accurately below, in Sec. IH-C. 



B. BOW-CDW transition 

In addition to proving the existence of the BOW phase, 
we have studied in detail the nature of the continu- 
ous BOW-CDW transition for two different values of U 
(U < Um)- For([/,y) = (C/c,Vc), i.e., on the BOW-CDW 
phase boundary, the real space staggered charge and ki- 
netic energy correlation functions fall off algebraically as 

{n,n,+r){-iy - r-\ 
{{K,K,+r) - {K,?){-ir - r-\ (10) 
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Based on conformal field theory calculations for similar 
phase transitions in ID spin systenis,E3 the exponent 77 
can be expected to depend on {Uc, K) but should be the 
same for both the CDW and BOW correlations. This 
gives the finite-size scaling of the structure factor and 
the susceptibility at the critical point: 



Xcdw,bow{tt) ^ N^^"^. 



(11) 



With a spin gap but no charge gap, as was demonstrated 
above, we expect tiie critical state to be of the Luther- 
Emery liquid type.cS The exponent 77 is then related to 
the Luttinger liquid parameter Kp hy rj — I — Kp. 

Fig. presents plots of ln[xc£)vi/] and ln[x_BOw] versus 
ln[iV] for [/ = 4 and three different values of V around 
the critical point, which as discussed above should be 
close to 2.16. The data points for y=2.16 indeed fall al- 
most on straight lines, indicating critical scaling for both 
the CDW and BOW fluctuations. The value of the criti- 
cal exponent ry, obtained from the slope of the V — 2.16 
curves for both xcdw and xbow is 77 w 0.5. The scaling 
of the structure factors. Scow and Sbow, aX V = 2.16 
is also consistent with 77 « 0.5. It is, however, difficult 
to extract a precise value for 77 from this finite-size scal- 
ing, due to subleading corrections to the scaling, as well 
as effects from the fact that the [/, V point studied is 
not exactly on the phase boundary. As was discussed in 
Sec. III-A, the Luttinger liquid parameter Kp can also be 
extracted from the g — > limit of ScDw{q) / q_- This is in 
general a more accurate method, since the convergence 
with system size is faster for the subleading 1/r^ contri- 
bution tp-.tlip correlation function which this estimator 
accesses Fig. || shows results for J7 = 4 and U — i 
and the respective critical V-values. The q behavior 
gives Kp = 0.44 ± 0.01 for [/ = 4, i.e., 77 = 0.56 ± 0.01, 
which hence is consistent with the finite-size scaling of 
the q = TT quantities. For J7 = 3 we,abtain Vc = 1.65, 
in agreement with Nakamura's result ,E3 and the critical 
exponent rj = 0.47 ± 0.01. 



C. SDW-BOW transition 

The SDW-BOW transition is marked by the opening of 
a spin-gap in the electronic energy spectrum. As argued 
by Nakamura, it is a quantum phase transition of the 
Kosterlitz-Thouless type and therefore the gap opens up 
exponentially slowly. This makes it difficult to determine 
the phase boundary numerically. The numerical data is 
affected by large finite-size effects that persist up to very 
large system sizes. As discussed in Sec. Ill A, the most 
reliable evidence of the existence of a spin-gap is obtained 
from the behavior of Ssdw{<i) I <1 as 5 ^ 0. In practice, 
an asymptotic value of ttSsdw{q)/q < 1 as g ^ in any 
(large) system is an indication of the presence of a spin- 
gap in the thermodynamic limit. This allows us to detect 
the presence of very small spin gaps. Fig. ^ shows the 
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FIG. 7. ln[xcDw{TT)] and Hxbow{tt)] vs ln[iV] for C/=4 
and different values of V near the critical point. The dashed 
lines are fits to the V — 2.16 data. 
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BOW-CDW boundary. 
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FIG. 9. SsDw{<l) /q vs. q for U — 4 and 4 difFerent values 
of V around the SDW-BOW boundary. 

behavior of S s dw {<i) 1 1 for [7 = 4 and difFerent valuta 
of V . In the gapless region, logarithmic correctionsE3 
makes it difficult to observe the-approach to 1 as q ^ 0. 
In analogy with spin systems J23 we expect the leading 
log corrections to vanish at the point where the spin gap 
opens, and therefore exactly at the critical point there 
should be a clear scaling to 1. An apparent reduction of 
the log correction is indeed seen in Fig. ^ as 1^ is increased 
towards « 1.88. Based on the results, we estimate the 
SDW-BOW boundary to be at V = 1.89 ±0.01 at [/ = 4. 
This is slightly higher than Nakamura's critical value V = 
1.82 for this U . We believe the difference is due to non- 
asymptotic finite-size effects in the exact diagonalization 
calculation, which used system sizes only up to = 14. 
Hence, we find that the BOW phase exists in a slightly 
smaller, while still significant, region of the phase space. 

D. First-order SDW-CDW transition 

For U > Um, the transition is a discontinuous (first- 
order) direct SDW-CDW transition with no intervening 
BOW phase. Fig. |l^ shows the V dependence of the 
CDW order parameter, the total energy, and the kinetic 
energy across the phase bouiid^ry for U = 8, which ac- 
cording to previous studiesB^lliJ'EJ should be well within 
the regime of first-order transitions. The characteristics 
of a first-order transition are indeed quite apparent. The 
order parameter and the kinetic energy change rapidly 
at the transition point, Vc ~ 4.14. The finite-size effects 
diminish with increasing N as the results approach the 
limiting behavior of a discontinuity in the order param- 
eter and the kinetic energy in the thermodynamic limit. 
The total energy remains continuous, but there is a clear 
break in slope at the transition. 

The size dependence of the BOW order-parameter is 
shown in Fig. |ll]. It becomes considerably smaller in- 
side the CDW phase than before the transition. This is 
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FIG. 10. The behavior of the CDW order parameter, 
the kinetic energy and the ground state energy across the 
SDW-CDW transition for various system sizes and U — 8. 

expected, since in the SDW phase, but not in the CDW 
phase, there should be power-law decaying BOW correla- 
tions. The BOW order parameter decays rapidly with the 
system size, however, confirming that there is no long- 
range BOW for this U > Um- 

The behavior with increasingly sharp discontinuities 
seen in Fig. |l^ and |l^ indicates a first-order transition 
due to an avoided level crossing. Note that with increas- 
ing chain length the CDW order parameter approaches 
its thermodynamic value from above for V < Vc and 
from below for V just above Vc- The curves for different 
system sizes cross one another in the neighborhood of 
V — Vc and then once again for a higher V . The second 
crossing point moves down towards the first one as iV in- 
creases, whereas the first crossing does not change much 
with V and appears to be a good criterion for locating 
the transition point. 

The two curve crossings can be understood as follows: 
In a transition caused by an avoided level crossing, a 
crossing of the order parameter curves close to the critical 
coupling (approaching the critical coupling as iV — > oo) 
can be expected since the low-energy levels correspond- 
ing to an ordered and disordered state swap characters 
within a a parameter range y± Ay (iV), with Ay (TV) 
as ^ oo. This behavior is seen clearly in Fig. |l^. The 
finite- Af ground state starts to develop CDW characteris- 
tics at Ay (TV) and thus, for a fixed V < Vc, the CDW 
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FIG. 11. The behavior of the BOW order parameter across 
the SDW-CDW transition for various system sizes and U — 8. 



order parameter decreases with increasing N. An anal- 
ogous argument for fixed V > Vc close to Vc gives that 
the in this case the CDW order must increase with in- 
creasing N. On the other hand, for V ^ Vc the real-space 
CDW correlations are enhanced at short distances (in the 
same way as the BOW correlations shown in Fig. ^) and 
for small system sizes there is also some enhancement of 
the long-distanjee correlations due to the periodic bound- 
ary conditions .cZl Hence, one can expect the CDW order 
parameter, when defined and measured in terms of its 
squared expectation Eq. (||), to again decrease with N 
for V ^ Vc and this explains the second crossing of the 
order parameter curves seen in Fig. 00. 



E. Multi-critical point 

Although the existence of the tricritical point (which, 
in vew of the existence of the BOW phase, we refer 
to as the multi-critical point) separating the first-order 
and continuous transition to the CDW state has long 
been known, its location in the (U, V) plane has not 
previously been dfttermined accurately using large sys- 
tem sizes. HirschB^EI estimated a value of Um—f= 3 using 
world line Monte Carlo. Cannon and Fradkinll2l obtained 
Um = 1.5 using field theory techniques and world lines. 
Later Cannon, Scalettar and FradkinEll obtained a value 
of Um = 3.5 — 5 using finite-size scaling of Lanczos re- 
sults. Using_a combination of bosonization and RG tech- 
niques, VoitO obtained Um — 4.76. However, as Voit also 
pointed out, the validity of bosonization and RG, which 
are applicable in the limit U,V — *■ 0, for intermediate 
values of the parameters is a priori questionable. 

By using larger system sizes and an alternative crite- 
rion to distinguish between a continuous transition and a 
first-order level crossing transition, we have obtained an 
estimate of the multi-critical point that we consider more 
accurate and reliable than the previous estimates. In con- 
trast to most previous numerical studies, our method is 
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FIG. 12. CDW order parameter vs V across the 
BOW-CDW boundary for several system sizes near the 
multi-critical point. The dashed line shows the position of 
Vc for the respective U. Statistical errors are smaller than 
the symbols. 
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not based on plotting histograms of the order parameter, 
although we will also present such histograms in the next 
section. In this section we first exploit the qualitatively 
different finite-size dependence of the growth of the or- 
der parameter close to the transition above and below 
the multi-critical point. 

For fixed U, the order parameter curves for different 
system sizes cross each other at or very close to the crit- 
ical point {V = Vc) in the case of a first-order transition, 
as discussed above in Sec. III-D. Such a crossing can- 
not occur at a continuous transition, where instead there 
should be finite-size scaling governed by Eq. (pi]). This 
qualitative difference in the finite-size dependence of the 
order parameter close to the transition point above and 
below the multi-critical point (Um,Vm) leads us to ex- 
pect that in the neighborhood of this point, curves of the 
order parameter for different chain lengths will closely 
coincide with one another close to V = Vc, and Um is the 
point at which the curves barely touch each other. When 
the system size becomes sufficiently large one can also di- 
rectly observe discontinuities developing when U > Um, 
in the order parameter as well as in other quantities, as 
in Fig. |o|. In practice this criterion, or any other crite- 
rion known to us, cannot be expected to be useful very 
close to the multi-critical point, where the transition is 
only weakly first-order and very large lattices are needed 
to detect discontinuities developing from avoided level 
crossings. 

Fig. na shows the finite-size dependence of the CDW 
order parameter across the transition for three different 
values of U. For U — 4.2, only the = 16 curve crosses 
the other curves, and this occurs far from the critical 
point (as determined using the peak in the charge stiff- 
ness, as discussed in Sec. III-A). The non-crossing for 
larger system sizes show that the transition must be con- 
tinuous at this U. For U = 5.2 all curves show crossing 
behavior and a discontinuity can also be seen develop- 
ing for the largest system size, i.e., the transition is here 
of first order. The curves for U — 4.6 closely follow 
the expected behavior at the multi-critical point, with 
the curves for the largest systems barely touching each 
other. Based on data also for other values of U we esti- 
mate the multi-critical point to be {Um = 4.7±0.1, Vm — 
2.51 ± 0.04). Xbis agrees very well with Voit's estimate 
{Um = 4.76).L3 However, it is not clear whether this 
agreement is fortuitous or whether there is some under- 
lying symmetry that renders bosonization and RG (that 
assume U,V ^ 1) applicable close to the multi-critical 
point. 

F. CDW order parameter histograms 

Previous studies of the multi-critical point have ex- 
ploited the existence of a 3-peak structure in the distri- 
bution of the CDW order parameter for a discontinuous 
SDW-CDW transition in the vicinity of the critical point 



and its absence at a continuous transition.B Outside the 
CDW phase, the distribution of the CDW order param- 
eter is peaked around zero. For a continuous transition 
to a CDW state this peak splits into two (corresponding 
to positive and negative values of the order parameter), 
which gradually move apart from each other inside the 
CDW phase. In a first-order transition, on the other 
hand, the order parameter takes a non-zero value imme- 
diately as the CDW phase is entered and hence the two 
peaks emerge already separated from each other. Fur- 
thermore, at the phase boundary the CDW phase coex- 
ists with the competing phase, and this is reflected as 
a central peak remaining in the CDW order parameter 
distribution. The position of the multi-critical point can 
then in principle be obtained by locating the point where 
the 3-peak structure first appears. In practice, the accu- 
racy of this method is limited by the fact that the dis- 
continuity is very small for a first-order transition close 
to the multi-critical point and very large system sizes are 
then needed to observe the three peaks. This problem 
is, of course, common to all methods for distinguishing 
between a continuous and weakly first-order transition. 

In his early QMC study, Hirsch observed a 3-peak 
structure even for U as small as 3 and therefore cori-i 
eluded that the transition is of first order already there.El 
For larger U, an unexplained 4-peak structure was seen. 
We have repeated histogram calculations for the lattice 
size = 32 studied by Hirsch. In Fig. |lj we show results 
for U = 6, V — 3.15, where a 4-peak structure was seen 
in the earlier calculation .El We only find a central peak, 
which show that the system is not in the CDW state 
for these parameters. There are, however, already signs 
of side peaks developing, which shows that the system 
is close to the CDW phase. The significant differences 
with the earlier result could partially be errors due to 
the Trotter decomposition used in the world-line simula- 
tion method. Temperature effects are only minor, as also 
shown in Fig. |l^. At /3 = 8, which was used in Ref. ||, the 
histogram is only slightly more sharply peaked than at 
/3 = 16 and 32. Most likely, the simulation giving the 4- 
peak structure was not sufficiently long, as it consisted of 
only lO'' Monte Carlo steps. cl Even with the more efficient 
SSE algorithm used in the present work, we find that the 
autocorrelation times are quite long close to the first or- 
der transition (see Appendix B) and short simulation can 
produce incorrect order parameter histograms similar to 
those shown in Ref. |8[ For the histograms shown here, 
on the order of 10'' — 10^ SSE Monte Carlo steps were 
used. 

In Fig. |l^ we also show results for several values of 
V across the phase transition. A clear 3-peak structure 
(i.e., three peaks in the range mcDW G l]i of which 
we only show the positive part) with peaks of almost the 
same heights can be seen for V = 3.165. In Fig. |lj we 
show results for N = 64. At U — 6, the 3-peak struc- 
ture appears for V ~ 3.156, i.e, at a value slightly lower 
than for the iV = 32 system. The size of the V region 
in which three peaks can be observed is also significantly 
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FIG. 13. CDW order parameter distributions for 32-site 
systems at [/ = 6. Upper panel: Dependence on the inverse 
temperature P at V = 3.15. Lower panel: Dependence on V 
around the first order phase transition. Statistical errors are 
of the order of the size of the symbols. 



smaller, reflecting the sharpening of the first order transi- 
tion cause by an avoided level crossing. At U = 5, which 
we have argued above should be close to but above the 
multi-critical point, we do not observe three peaks. The 
histogram becomes very flat for an extended range of 
mcDW, however, and the side peak emerges at a finite 
value of mcDW- This is consistent with the transition 
still being of first order a,t V = 5. Going to still lower 
^-values, the peak just becomes narrower, and it is not 
possible to definitely conclude this way when the transi- 
tion becomes continuous. 



IV. DOPED SYSTEMS 

An interesting question naturally arises from the ex- 
istence of the Nakamura BOW phase: Can the EHM 
model support a soliton lattice when doped slightly 




FIG. 14. CDW order parameter histograms for A'^ = 64 
systems close to the phase transition. 



away from half-filling? Such a state exists in the Su- 
Schrieffer-Heeger model ii*-tli©iadiabatic limit (with clas- 
sical, "frozen" phonons),EZrL3 and also wlien—electron- 
electron interactions are taken into account .ljO In these 
models, the quantum nature of the phonon field is not 
taken into account, however. It is known that the dimer- 
ized state at half-filling survives even in the presence 
of quantum fluctuations, j-at least up to a critical value 
of the phonon frequency.E3 However, to our knowledge, 
there have been no reliable numerical calculations ad- 
dressing the stability of the soliton lattice in the pres- 
ence of fully quantum mechanical phonons. The Naka- 
mura BOW is similar to a dimerized lattice with quantum 
fluctuations, and hence a study of its evolution with hole 
doping can give insights also into the quantum phonon 
problem. There are_also unresolved issues regarding the 
doped CDW state.Ej In order to investigate the evolu- 
tion of the long-range ordered states upon doping, we 
have studied the EHM model also away from half-flUing, 
focusing on two parameter values in which the half-filled 
system is in the BOW (using [/ = 4, F = 2.14) or CDW 
phase (using U = A, V = 2.5). 

We first discuss the effects of doping on the spin and 
charge gaps in the half-filled CDW and BOW ground 
states. As in previous sections, we make use of the be- 
havior of the static structure factor in the limit g ^ 0. 
Figs. ^ and |l^ show ttSc dw (m) /q and nSs dw (?) /g as a 
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FIG. 15. The static charge (upper panel) and spin (lower 
panel) structure factor divided by the wave-number q as a 
function of q for 256-site chains at different doping fractions 5. 
For these parameter values ()7 = 4, 1/ = 2.14), the half-filled 
system (5 = 0) is in the BOW state. 
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FIG. 17. Static CDW and SDW susceptibilities at different 
doping levels for U — 4, V = 2.5 for a A'' = 256 chain. The 
inset shows the xcDw{2kp) peaks on a more detailed scale. 



function of q for a range of doping levels, both for param- 
eters where the half-filled system is in the BOW phase 
(Fig. |l5|) and in the CDW phase (Fig. |l|). From the data 
we conclude that upon doping away from half- filling, the 
charge gap vanishes immediately whereas the spin gap 
survives. This is true for both the CDW and BOW parent 
stateSj-This behavior is characteristic of a Luther-Emery 
liquidjlj in which the charge sector can be described in 
terms of a Luttinger liquid and the spin sector is gapped. 
The limiting value of nSc (g) as g — > indicates that 
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FIG. 16. Same as Fig. ^ior U = 4,V = 2.5, where the 
half-filled system is in the CDW phase. 



the Luttinger liquid exponent Kp « 0.5 in both the cases, 
with only a weak dependence on the doping level for the 
parameters considered here. Note the cross-over behav- 
ior occurring in the charge structure at g ~ 2ttS = Akp 
in Fig. 1^ (which is accompanied by a peak in the-ror- 
responding susceptibility, as will be shown below) Jld re- 
flecting a weak repulsion between dopant holes. No cross- 
over in the charge structure is seen in Fig. where the 
parent state is a BOW. 

Fig. ^ shows the variation of the ground state static 
susceptibilities for several doping levels in a chain of 
length N — 128 for the parameters U = 4:,V = 2.5. For 
6 > the charge susceptibility converges to a non-zero 
value as q 0, again showing the absence of a charge 
gap. Very strong 2kF peaks are evident, and weaker Akp 
peaks are also clearly visible. The 2kF peaks diverge with 
the system size whereas the Akp peaks are non-divergent, 
in accord with the Luther-Emery picture. For a Luther- 
Emery liquid, the charge correlations decay with distance 
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FIG. 18. Finite-size scaling of the static charge susceptibil- 
ity at g = 2kF for a systems with ?7 = 4,l/ = 2.5ata doping 
level of 6.25%. A slope of 1.5 is shown by the dashed line. 

r as r~-^p,Q which gives XcDwi'^kp) N'^~^p for the 
finite-size scaling of the corresponding 2kp susceptibility. 
Fig. |l^ shows the size dependence for 6 — 0.0625 on a 
log-log scale. For system sizes N > 64, the data is seen to 
fall on a line with slope w 1.5, consistent with the value 
Kp fa 0.5 extracted above. 

As a further test of the Luthcr-Emcry liquid nature of 
the ground state away from half-filling, we have studied 
the real-space charge and bond correlations as a function 
of distance. Fig. |l^ shows the charge correlation for two 
different system sizes at a dopant concentration of 6.25% 
and interaction parameters U = A and V — 2.5. The 
ground state at half-filling is a CDW, and for the doped 
system we find solitonic features with alternating A and 
B phases separated by domain walls or kinks. The corre- 
lation decay with distance, however, and there is no real 
soliton lattice. In fact, the decay of the magnitude of the 
peaks is well approximated by an envelope curve of the 



form y - 



-0.5 



and hence also these data are consistent 



with a Luther-Emery state with Kp « 0.5. 

Fig. ^ shows a similar plot of the real-space bond- 
order correlation for U = A and V = 2.14 at the same 
dopant concentration and for the same system sizes. The 
ground state of the half-filled system for this choice of 
parameters is here a BOW, and away from half-filling, the 
dominant correlation are still of bond-order type. Once 
again, the ground state of the doped system has solitonic 
features with an algebraic decay of the magnitude of the 
peaks. As in the previous case, the decay is consistent 
with a Luther-Emery liquid with Kp w 0.5. 



V. SUMMARY 

To summarize, we have studied the ID EHM using the 
SSE method incorporating an efficient operator-loop up- 
date and a "quantum parallel tempering" scheme. Our 



0.2 



0.1 



0.0 



-0.1 



-0.2 



0.2 



0.1 



0.0 



-0.1 



-0.2 




10 



30 



\ 



50 



N=256 



60 




80 100 



FIG. 19. Real-space charge correlations vs distance at a 
doping level 5 = 0.0625 fo system sizes = 128 and 256 at 
U = A,V = 2.5. 
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FIG. 20. Bond-order correlations at 5 = 0.0625, 
f/ = 4, V = 2.14, for system sizes = 128 and 256. 
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results confirm the surprising prediction0 of the exis- 
tence of a novel long-range-ordered BOW phase between 
the well-known CDW and SDW phases in the ground 
state phase diagram for small to intermediate values of 
the on-site interaction U {U < Um)- We have presented 
several ways to detect the spin and charge gaps expected 
in the BOW phase and have also probed directly the 
BOW correlations and concluded that true long-range 
order develops. We have studied a few points on the 
BOW-CDW phase boundary and obtained a very goad, 
agreement with Nakamura's level crossing predictionE£l 
for the location of this phase boundary. For the SDW- 
BOW phase boundary, our results indicate a higher criti- 
cal V for fixed U than given by the level crossing method 
and thus over-all a slightly smaller size of the BOW 
phase. Our results are for significantly larger systems 
than in the previous study and it is not surprising that 
the finite-size effects in the level crossings can be large for 
the SDW-BOW transition since the spin gap opens ex- 
ponentially slowly in this Kosterlitz-Thouless transition. 
An over-estimation of the size of the BOW phase from 
the level crossings is also apparent considering that our 
estimated multi-critical point is well within the BOW 
phase of Nakamura's phase digram. Since our BOW- 
CDW phase boundaries agree, this indicates problems 
with the scaling of the exact SDW-BOW level crossings 
close to thej_multi-critical point, as was also mentioned by 
Nakamura.Ej For large values of [/ ([/ > Um) the tran- 
sition is discontinuous (first-order) . We have shown that 
curves of the CDW order parameter across this bound- 
ary for different system sizes cross each other twice, and 
explained this behavior in terms of an avoided level cross- 
ing. We have also used the curve crossings as a means to 
locate the position of the multi-critical point with greater 
accuracy than previously attained. Our estimate for the 
multi-critical point is [/,„ = 4.7 ± 0.1, Vm = 2.51 ± 0.04. 

We have also studied systems doped slightly away from 
half-fiUing. We find that both the doped CDW and BOW 
states give rise to ground states of the Luther-Emery 
type, i.e., the quantum fluctuations do not allow the for- 
mation of true soliton lattices. Based on the fact that the 
BOW state is very similar to the dimerized ground state 
of models with finite- frequency (non-adiabatic) phonons, 
we conjecture that the soliton lattice is also unstable to 
arbitrarily weak quantum fiuctuations in these models, 
unless two- or three-dimensional couplings are taken into 
account. 

After our completion of the nuinerical calculations at 
half-filling Tsuchiizu and Furusakic3 presented a weak- 
coupling g-ology calculation taking into account second- 
order corrections to the coupling constants. They ob- 
tained a phase diagram in very good quantitative agree- 
ment with ours, including the location of the multi- 
critical point. 
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APPENDIX A: OPERATOR-LOOP UPDATES IN 
THE SSE METHOD. 

The SSE approach has been discussed in several 

papersJl3il3. We here start with a brief pview as a basis 
for introducing the operator-loop updatetJ in the context 
of fermion models. 

To implement the SSE method, the Hamiltonian (^ is 
written, upto an additive constant, in the form 



N 



(Al) 



6=1 



where b is the bond connecting the sites b and b+1, N is 
the length of the chain, and the operators Ha.b, a — 1,2,3 
are defined as 



iJi.b = C-|(nb,T-i)(7ib,i-5) 
- f ("6+1,1 - ^ ' 



H2,b 
H3,b 



-1,T - 2A"6+l,i 

V{nb - l)(nb+i - 1), 
^(4+i,iC6a +h.c.), 
<(4+i,t'^''-T +h.c.). 



(A2) 



The constant C shifts the zero of the energy and is 
chosen to ensure a non-negative expectation value for 
Hi^i, (needed in order to ensure a positive definite ex- 
pansion of the partition function). Introducing a basis 
{|a)} = {lCi,C2,---,Civ)}, where e {0, t, i, Ti} denotes 
the electron state at the site i, the partition function 
Z=Tr{e~^^} can be expanded in a Taylor series as 



EEE 

a n=0 S„ 



/3" 



(A3) 



where Sn denotes a sequence of index pairs defining the 
operator string np=i Ha^^b^: 



Sn = [a,b]i[a,b]2 ■ . . [a, 6], 



(A4) 



where we use the notation [a, b]p — [up, bp] and a £ 
{1, 2, 3}, 6 e {1, ... , N}. In order to construct an effi- 
cient updating scheme, the Taylor series is truncated at 
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a self-consistently determined power L, large enough to 
cause only an exponentially small, completely negligible 
error (L ~ l3\E\, where E is the total internal energy; 
for details see Refs. [TtUT^ ). We can then define a sam- 
pling space where the length of the sequences is fixed, 
by inserting L ~ n unit operators, denoted by Hq q, into 
each sequence. The terms in the partition function must 
be divided by (;^) in order to compensate for the differ- 
ent ways of inserting the unit operators. The summation 
over n is then implicitly included in the summation over 
all sequences of length L. The partition function takes 
the form 



a Sl 



/?"(L-n)! 
LI 



(A5) 



where the operator-index pairs [a, b]p now have a £ 
{1,2,3} and be {1, . . . ,iV} or [a,b]p = [0,0]. For con- 
venience, we introduce a notation for states obtained by 
the action of the first p elements of the operator string 
Sl- 



(A6) 



For a nonzero contribution to the partition function, 
\a{L)) = |a(0)). 

A Monte Carlo scheme is used to sample the configu- 
rations (a, Sl) according to their relative contributions 
(weights) to Z. i-Tte sampling scheme consists of two 
types of updates Jljilj referred to as diagonal update and 
operator- loop updates. The diagonal update involves lo- 
cal substitutions of the form [0,0]p ^ and is at- 
tempted consecutively for every p € {1,...,L} in the 
sequence for which [a, b]p = [0, 0]p or [1, b]p. The updates 
are accepted with probabilities 



P{[0,0]p 



P{[l,b]p^[0,0]p) 



NfiMiM 

L — n 
L-n + 1 

NfiMiM ' 



(A7) 



where 

Ma,b{p) = {Cb{p)Xb+l{p)\Ha.b\Cb{p-l)Xb+l{p-l)) 

(A8) 

is a matrix element on bond 6, which in this case is diag- 
onal (a — 1). Only a single state \a{p)) is stored in the 
computer during the diagonal update. When off-diagonal 
operators are encountered during the successive scanning 
of the operator string, the corresponding electron states 
are updated so that the information needed for evalua- 
tion of the probabilities (A7) is always available when 
needed. 

The operator-loop update has been discussed in de- 
tail in Rcf. [13 in the context of spins. Here we present 



the construction of loops for fermio ns. As explained in 
Ref. the matrix element in Eq. ( [A5| ) can be graph- 
ically represented by a set of n vertices (corresponding 
to the n non-unit operators in Sl) connected to one an- 
other by the propagated electron states. Each vertex has 
four "legs" with electron states \C,i{p—l),Q+i[p—l)) and 
\Ci{p) : Ci+i{p)) before and after the action of the associ- 
ated Hamiltonian operator Hap^bp- There are 32 allowed 
vertices - 16 diagonal ones and 8 each associated with 
the off-diagonal H2.b and H^ b [see Fig. ^(a)]. A config- 
uration (ujSl) is completely specified by the leg states 
of the n vertices — except for sites that do not have any 
operators acting on them. 

To carry out the operator-loop update, the linked list 
of the n vertices is first constructed. In addition to the 
electron states at the legs of each vertex, the list also con- 
tains the addresses (i.e., the location in Sl) of the next 
vertex and the corresponding leg that each leg is con- 
nected to. The loop construction begins with randomly 
choosing a vertex and an "entry" leg. The electron state 
at the entry leg is changed to one of the 3 other allowed 
states chosen at random. Next an "exit" leg is chosen 
(following a procedure described below) and its associ- 
ated electron state is updated so that the new leg states 
constitute an allowed vertex [see Fig. |2^(b)]. The exit 
leg will be linked to a leg of another vertex (or, if there is 
only one operator in the configuration which acts on the 
site in question, another leg on the same vertex) and this 
will be the entry leg for the next vertex. The electron 
state at this new entry leg is then updated to match the 
state at the exit leg of the previous vertex. A new exit 
leg is then chosen following the same procedure. This is 
repeated until the exit leg from a vertex points to the 
starting point of the loop, which implies that the loop is 
closed and a new allowed configuration has been gener- 
ated. 

To choose an exit leg — given a vertex, an entry leg 
and the updated electron state at the entry leg — all 
the legs can be considered in turn and attempts made 
to update the associated electron state so that the new 
leg states constitute an allowed vertex. Because of spin 
an charge conservation on the vertices, at a given exit 
leg there is at most one possible update of the electron 
state that can lead to an allowed vertex. Hence, the exit 
leg uniquely determines the new vertex and the proba- 
bility of choosing a given leg should be proportional to 
the weight of the new vertex, i.e., a matrix element of 
the form (A?), which in this case can be either diagonal 
or off-diagonal. In practice, a fast selection of an exit leg 
and updating of the vertex state is achieved using two 
pre-generated tables. The first one contains the cumu- 
lative probabilities of the 4 exit legs given an entrance 
leg, the old vertex state, and the new state at the en- 
trance. The second table contains the new vertex states 
corresponding to the updated entrance and exit legs. 

A special case occurs if the initial update at the entry 
leg of the first vertex of a loop is a spin-flip, i.e. the elec- 
tron state changes from "f to | or vice versa. In this case. 
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FIG. 21. (a) A few allowed vertices. The solid lines denote 
the diagonal Hamiltonian operator, the dashed and dotted 
lines denote the hopping operators for the up and down spins 
respectively. The lower legs denote the states ^i{p — 1) and 
Ci+i(p— 1) while the upper legs denote Ci(p) and (b) 
An example of a vertex update. The entrance leg is the lower 
left leg of the vertex, as indicated by the dot. The electron 
state at the entrance leg, |, is changed to || in this particu- 
lar update. Given that, the 3 possible resulting vertices are 
shown. The corresponding exit legs are denoted by open cir- 
cles. Exit at the upper right leg does not result in an allowed 
vertex in this case. 



the vertex weight does not change when updated and as 
a consequence the "bounce process" , where the exit leg is 
the same as the entrance leg, does not have to be included 
in the loop construction. The loop then becomes deter- 
mirusjtic, i.e., there is a unique exit leg given the entrance 
leg.ll3 This is similar to the "loop-exchange" algadthm 
proposed in the context of the world-line method.c3 

A full Monte Carlo updating cycle (MC step) consists 
of a diagonal update, followed by the construction of a 
linked vertex list. Next a number of operator- loop up- 
dates are carried out and finally the vertices are mapped 
back into a corresponding sequence Sl- The loop up- 
date typically also implies changes in the stored state 
|a) = |a(0)), as some of the vertex legs (links) span across 
the periodic boundary in the propagation direction. The 
number of up and down electrons can be changed by the 
operator-loop update, as can the spatial winding num- 
bers, and the algorithm is hence fully grand canonical. 
Note that at high and moderately low temperatures there 
are typically some sites of the system which has no ver- 
tices associated with them. The states on these sites can 
be randomly changed, since they have no affect the con- 
figuration weight. 

The number of loops constructed for every MC step is 
determined such that on an average a total of ^ L ver- 
tices (we typically use 2L) are visited. The truncation 
L and the number of loops are adjusted during he equi- 
libration part of the simulation and are thereafter held 
fixed. L is determined by requiring that the highest n 
reached during equihbration is at most 70 — 80% of L. 

In certain parameter regions the length of a loop can 
sometimes become extremely long before it closes — in 
practice, it may even never close. It is therefore necessary 
to impose a maximum length, beyond which the loop con- 



struction is terminated and a new starting point is chosen 
(typically, we use « SOL for this cut-off length). In order 
to reduce the likelihood of the next loop also exceeding 
the termination length it has proven useful to carry out a 
diagonal update before starting the next loop. The loop 
termination does not violate detailed balance and does 
not cause any systematical errors in the results. In most 
cases, incomplete loop termination occurs so infrequently 
that it does not adversely affect the simulation. In anal- 
ogy with Ref. |4j, where a scheme (there called "worm" 
update) similar to the operator-loops considered here was 
first introduced within the continuous worl-line represen- 
tation, the end points of the loop during construction 
can be related to the single-particle Green's function of 
the system and hence the tendency for loops to become 
exceedingly long for some parameter values must be re- 
lated to some physical properties of the system. This 
issue should be studied further. 

Estimators for the various structure factors andj-wis^ 
ceptibilities have been discussed in previous articles O'tS 
Here we only note that the charge and spin stiffness con- 
stants, Eq. (|^), can be expressed in terms of spin and 
charge current operators in analogy with the spin stiffness 
of the Heisenberg antiferromagnet previously discussed in 
Ref. |l^, leading to 



[{n 



Nf3 



(A9) 



where ' 



'■R,L 



are the number of kinetic energy operators in 



the SSE term propagating spin-cr particles in the "right" 
and "left" direction on the ring. Because of spin and 
charge conservation, the topological winding numbers 
{nfi — n1)/N can take only integer values. 

Although the operator-loop algorithm very signifi- 
cantly speeds up SSE simulations, in many cases reducing 
the autocorrelation function by orders of magnitude, the 
dynamics is still very slow in some parameter regions. For 
the extended Hubbard model studied here, problems with 
very long autocorrelation times occur in the long-range 
ordered BOW and CDW phases. The problems are par- 
ticularly severe for large systems close to the BOW-CDW 
phase boundary, where "trapping" in the wrong phase of- 
ten occurs. The slow dynamics in the BOW phase is 
illustrated in Figure |2^, which shows the simulation time 
dependence of SBOwiT^) and XBOwi^^) during a simula- 
tion of a 256-site system at /3 = 512 [Sbow has converged 
at this P but xbow is about 20% larger still at /? = 1024]. 
It is evident that the BOW autocorrelation time here is 
tens of thousands of MC steps. The BOW susceptibility 
exhibits a behavior where it sometimes takes very small 
values (less that 10~^ of the average value), but corre- 
sponding large fluctuations upwards do not occur, i.e., 
the distribution of the XBOwi''^) estimator for individ- 
ual configurations is very skewed. The structure factor 
exhibits a more symmetric distribution. This behavior 
can be understood as a consequence of the BOW ground 
state for a finite system being a symmetric combination 
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FIG. 22. The BOW structure factor and susceptibility for 
a 256-site system with (7 = 4 and V = 2.14 at inverse tem- 
perature j3 — 512. Results of six independent simulations are 
shown. Each point represents an average over a bin consisting 
of 10* Monte Carlo steps. 



of the two possible real-space symmetry-broken states. 
The symmetry is not broken in a finite system and the 
simulation is also not trapped in one of the real space- 
states. Hence, the wave function that is sampled in the 
simulation contains both the real-space states and the 
behavior seen in Fig. ^ indicates that individual con- 
figurations also contain both components, in such a way 
that transitions ( "tunneling" ) between the two real-space 
states can occur daring the SSE propagation (which can 
be simply relatedeS to a propagation in imaginary time), 
at least for some configurations. Tunneling can be in- 
ferred from the qualitatively different evolutions of the 
structure factor and the susceptibility in MC time. The 
susceptibility is an integral of the bond-order correlation, 
as in Eq. (jj), which in configurations where tunneling 
occurs can be much smaller than in configurations with 
no tunneling, since correlations between states with the 
same real-space configuration contributes positively but 
correlations between different states give a negative con- 
tribution. The structure factor, on the other hand, is an 
equal time correlation function and would not be much 
reduced by tunneling if the tunneling times are short. 
This explains the qualitatively different distributions of 
the xbow and Sbow measurements in Fig. ^ Evi- 
dently, the updating process is very slow in adding and 
removing tunneling events in the configurations, which 
maybe is not that surprising considering that the tunnel- 



ing is between two states with a discrete broken symme- 
try. These problems do not occur in SSE simulations of 
systems with a broken continuous symmetry, such as the 
two-dimensional Heisenberg model. 

The trapping and tunneling problems can be signif- 
icantly reduced by using tke parallel tempering scheme 
(or exchange Monte Carlo) ,E3~E3 which is discussed below 
in Appendix B. 



APPENDIX B: QUANTUM PARALLEL 
TEMPERING 

The "quantum parallel tempering" scheme is a 
straight-focwpjd generalization of the thermal parallel 
temperin^3~E3 method commonly used to equilibrate 
classical spin glass simulations. Our implementation 
amounts to running several simulations simultaneously 
on a parallel computer, using a a fixed value of U and 
different closely spaced values of V. Along with the usual 
Monte Carlo updates, we attempt to swap the configura- 
tions for processes with adjacent values of V at regular 
intervals, typically after every Monte Carlo step, accord- 
ing to a scheme that maintains detailed balance in the 
extended ensemble of parallel simulations. The probabil- 
ity of swapping the 1^- values of runs i and i + 1, which 
are running at Vi and l^+i, respectively, before the swap, 
is 

p f,r ^ . M Wi(Vi+l)Wi+l(Vi) ^ 

Pswapd^. l/.+i) = mm[l, ^^(v,)w.^,(Vi^r) ^' ^""'^ 

where Wi {V) is the weight of the ith simulation configu- 
ration evaluated with the coupling V. The swap proba- 
bilities for fixed AV = Vi+i—Vi decreases with increasing 
system size and decreasing temperature and hence AV 
and the range of F-values (if the number of processes is 
fixed) must be chosen smaller for larger system sizes. 

The computational effort required for the swapping 
process is very minor compared to the actual quantum 
Monte Carlo simulations. It is therefore useful to carry 
out several swap attempts of all pairs of neighboring sim- 
ulations between every MC step. Histograms containing 
the number of times each of the current configurations 
has "occupied" each T^-bin can then be constructed and 
used for adding the contributions of each configuration 
to all the y-bins. This can contribute to reducing the 
statistical error of measured quantities. 

To illustrate the advantage of quantum parallel tem- 
pering, we show two sets of data — obtained with and 
without the use of tempering — for a system undergoing 
a first-order transition. Fig. ^ shows the CDW order pa- 
rameter across the first-order SDW-CDW phase bound- 
ary at [/ = 8. The upper panel shows the data obtained 
from individual runs; the lower panel shows data for the 
same parameters obtained using tempering. The length 
of the individual simulations was 10^ MC steps for all V 
values, and this was also the number of steps performed 
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FIG. 23. The CDW order parameter across the SDW-CDW 
phase boundary for i7 = 8 (A*' = 64, /3 — 64). The upper 
panel shows data from individual runs. The lower panel shows 
the same data obtained using quantum parallel tempering, 
with two independent runs as indicated by the open and solid 
circles. 



by each process in the tempering runs. The improvement 
in the quahty of the tempering data is evident, especially 
close to the transition point where two of the individual 
simulations have relaxed into the wrong phases. The sta- 
tistical errors are hence severely underestimated due to 
the failure to equilibrate properly within the simulation 
time. The tempering error bars are also large at the phase 
transition, but in contrast to those of the individual sim- 
ulations they are accurate error estimates. The errors 
rapidly become much smaller as one moves away from 
the transition point. The effects of tempering are also 
favorable further inside the CDW phase, where several 
of the individual simulations are apparently affected by 
trapping in configurations with defects, where the order 
is reduced. 

The tempering acceptance rate during the run span- 
ning across the phase transition in Fig. 23 is shown in 
Fig. 
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There is a sharp reduction in the acceptance 
rate at the transition. This reflects the rapid change in 
the SSE configurations across the phase boundary, which 
implies that the configuration weights evaluated with V 
values from the "wrong" phase are likely to decrease and 



the swap according to the probability (Bl) to be rejected. 



Finally, we note that tempering, in general, is an ap- 
plication where a superlinear speed-up can be achieved in 
practice on parallel computers. In addition to doubling 
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FIG. 24. The tempering acceptance rate during the simu- 
lation across the first-order SDW-CDW phase boundary. 



the density of data points when the number of processes 
is doubled, the statistical errors are also reduced. Some- 
times the error reduction can be dramatic, but even in 
cases where there are no real problems with the dynam- 
ics of individual simulations the effects of tempering are 
often very favorable. 
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